home *** CD-ROM | disk | FTP | other *** search
/ SPACE 1 / SPACE - Library 1 - Volume 1.iso / program / 168 / library / math.c < prev    next >
Text File  |  1988-04-30  |  6KB  |  345 lines
  1. /* C library for MJC 2.0 - floating point and long arithmetic */
  2.  
  3. /* 
  4.  * my floating point routines 
  5.  *
  6.  * 32 bit floating point format
  7.  *    bit    function
  8.  *    31    sign
  9.  *    24-30    exponent (excess 64)
  10.  *    0-23    fraction
  11.  */
  12.  
  13. #define NULL    0L
  14.  
  15. #define EXCESS    64
  16. #define DIGITS    24
  17. #define MAXFRAC    0xFFFFFF
  18. #define MAXEXP    0x7F
  19. #define MAXBITS    32
  20.  
  21. #define sign(x) ((x >> 31) & 1)
  22. #define exp(x) (((x >> 24) & MAXEXP) - EXCESS)
  23. #define frac(x) (x & MAXFRAC)
  24. #define negate(x) (x ^ 0x80000000L)
  25.  
  26. /* long to float */ 
  27.  
  28. long
  29. _ltof(n) long n; {
  30.     int s;
  31.     long _fnorm();
  32.     if (s = (n < 0L)) n = -n;
  33.     return _fnorm(s, DIGITS, n);
  34. }
  35.  
  36. /* float to long */
  37.  
  38. long
  39. _ftol(u) long u; {
  40.     int e, s;
  41.     long f;
  42.     s = sign(u);
  43.     e = exp(u);
  44.     f = frac(u);
  45.     if (e <= 0) { /* no fractional part */
  46.         return 0L;
  47.     }
  48.     else if (e < DIGITS) {
  49.         f = f >> (DIGITS - e);
  50.     }
  51.     else if (e == DIGITS) {
  52.         /* already there! */
  53.     }
  54.     else if (e < MAXBITS) {
  55.         f = f << (e - DIGITS);
  56.     }
  57.     else    {
  58.         _ferr("ftol overflow\n\r");
  59.         f = 0x7FFFFFFF;
  60.     }
  61.     return s ? -f : f;
  62. }
  63.  
  64. long
  65. _fcmp(x, y) long x, y; {
  66.     return (sign(x) && sign(y)) ? y - x : x - y;
  67. }
  68.  
  69. long
  70. _fneg(x) long x; {
  71.     return negate(x);
  72. }
  73.  
  74. long
  75. _fsub(u, v) long u, v; { 
  76.     long _fadd();
  77.     return _fadd(u, negate(v)); 
  78. }
  79.  
  80. long
  81. _fadd(u, v) long u, v; {
  82.     int t, sw, eu, ev, ew;
  83.     long fu, fv, fw, _fnorm();
  84.     if (u == 0L) return v;
  85.     if (v == 0L) return u;
  86.     eu = exp(u);
  87.     ev = exp(v);
  88.     fu = sign(u) ? -frac(u) : frac(u);
  89.     fv = sign(v) ? -frac(v) : frac(v);
  90.     t = eu - ev;
  91.     if (t <= -DIGITS) { /* u << v */
  92.         fw = fv;
  93.         ew = ev;
  94.     }
  95.     else if (t <= 0) { /* u <= v */
  96.         if (t) fu >>= -t;
  97.         fw = fv + fu;
  98.         ew = ev;
  99.     }
  100.     else if (t < DIGITS) { /* u > v */
  101.         fv >>= t;
  102.         fw = fv + fu;
  103.         ew = eu;
  104.     }
  105.     else    { /* u >> v */
  106.         fw = fu;
  107.         ew = eu;
  108.     }
  109.     if (sw = (fw < 0))
  110.         fw = -fw;
  111.     return _fnorm(sw, ew, fw);
  112. }
  113.  
  114. long
  115. _fmul(u, v) long u, v; {
  116.     int sw, ew, eu, ev;
  117.     long fu, fv, fw, _fnorm();
  118.     if ((fu = frac(u)) == 0)
  119.         return 0L;
  120.     if ((fv = frac(v)) == 0)
  121.         return 0L;
  122.     eu = exp(u);
  123.     ev = exp(v);
  124.     sw = sign(u) ^ sign(v);
  125.     ew = eu + ev - 8;
  126.     fw = (fu >> 8) * (fv >> 8);
  127.     return _fnorm(sw, ew, fw);
  128. }
  129.  
  130. long
  131. _fdiv(u, v) long u, v; {
  132.     int sw, ew;
  133.     long _fnorm();
  134.     long b, fu, fv, fw;
  135.     if ((fu = frac(u)) == 0)
  136.         return 0L;
  137.     if ((fv = frac(v)) == 0)
  138.         return 0L;
  139.     sw = sign(u) ^ sign(v);
  140.     ew = exp(u) - exp(v) + 1;
  141.     fw = 0L;
  142.     b = 1L << (DIGITS - 1);
  143.     while (b) {
  144.         if (fu >= fv) {
  145.             fw += b;
  146.             fu -= fv;
  147.         }
  148.         b >>= 1;
  149.         fv >>= 1;
  150.     }
  151.     return _fnorm(sw, ew, fw);
  152. }
  153.  
  154. /* ascii to float, needs no other support routines */
  155.  
  156. long 
  157. atof(s) char *s; {
  158.     int e, i, c, sgn;
  159.     unsigned long ipart, fpart, f, bit, div;
  160.     long _fnorm();
  161.     /* get to the start of the digits */
  162.     while (*s && *s <= ' ')
  163.         s++;
  164.     if ((sgn = (*s == '-')) || *s == '+')
  165.         s++;
  166.     /* get the integer and fraction parts */
  167.     ipart = 0L;
  168.     while (*s >= '0' && *s <= '9')
  169.         ipart = ipart * 10L + (*s++ - '0');
  170.     fpart = 0L;
  171.     if (*s == '.') {
  172.         s++;
  173.         /* 9 digits of precision please */
  174.         for (i = 0; i < 9; i++) {
  175.             if (*s >= '0' && *s <= '9')
  176.                 c = *s++ - '0';
  177.             else    c = 0;
  178.             fpart = fpart * 10L + c;
  179.         }
  180.         div = 500000000;
  181.     }
  182.     /* normalize the integer part, then work in the fraction */
  183.     if (ipart) {
  184.         f = _fnorm(sgn, DIGITS, ipart);
  185.         e = exp(f);
  186.     }
  187.     else    e = f = 0L;
  188.     if (fpart && e < DIGITS) { /* room for the fraction part */
  189.         f = frac(f);
  190.         bit = 1L << (DIGITS - 1 - e);
  191.         /* add the fraction's bits to f */
  192.         while (bit && fpart) {
  193.             if (fpart >= div) {
  194.                 fpart -= div;
  195.                 f |= bit;
  196.             }
  197.             bit >>= 1;
  198.             div >>= 1;
  199.         }
  200.         if (fpart) f++; /* round up */
  201.         f = _fnorm(sgn, e, f);
  202.     }
  203.     /* now handle an exponent factor, if any */
  204.     if (f && (*s == 'e' || *s == 'E')) {
  205.         i = 0;
  206.         s++;
  207.         if ((c = (*s == '-')) || *s == '+')
  208.             s++;
  209.         for (i = 0; *s >= '0' && *s <= '9'; s++)
  210.             i = i * 10 + (*s - '0');
  211.         if (c) i = -i;
  212.         while (i > 0) {
  213.             e = exp(f);
  214.             f = frac(f);
  215.             f = _fnorm(sgn, e, f * 10L);
  216.             i--;
  217.         }
  218.         while (i < 0) {
  219.             e = exp(f);
  220.             f = frac(f);
  221.             f = _fnorm(sgn, e, f / 10L);
  222.             i++;
  223.         }
  224.     }
  225.     return f;
  226. }
  227.  
  228. long
  229. _fnorm(s, e, f) int s, e; unsigned long f; {
  230.     if (f == 0) {
  231.         return 0L;
  232.     }
  233.     else    {
  234.         while (f & 0xFF000000L) {
  235.             f = (f >> 1) + (f & 1);
  236.             e++;
  237.         }
  238.         while ((f & 0x800000L) == 0) {
  239.             f = f << 1;
  240.             e--;
  241.         }
  242.     }
  243.     if (e >= EXCESS) {
  244.         _ferr("float overflow\n\r");
  245.         return ((long)s << 31) | 0x7FFFFFFFL;
  246.     }
  247.     else if (e < -EXCESS) {
  248.         _ferr("float underflow\n\r");
  249.         return 0L;
  250.     }
  251.     /* pack the result */
  252.     return     ((long)s << 31) | 
  253.         (((long)(e + EXCESS) & MAXEXP) << DIGITS) | 
  254.         (f & MAXFRAC);
  255. }
  256.  
  257. _ferr(s) char *s; { while (*s) trap(1, 2, *s++); }
  258.  
  259. /* 
  260.  * stuff for long binary ops (both signed and unsigned)
  261.  *    a and b are the left and right arguments to the binary operator
  262.  *
  263.  * div and mod are done by long division, shift b up until >= a, then
  264.  * back down, subtracting when appropriate
  265.  *
  266.  * mul is done by shifts and adds
  267.  */
  268.  
  269. long
  270. _lmul(a, b) long a, b; {
  271.     int neg = 0;
  272.     long _ulmul();
  273.     if (a < 0L) { neg++; a = -a; }
  274.     if (b < 0L) { neg++; b = -b; }
  275.     a = _ulmul(a, b);
  276.     return neg & 1 ? -a : a;
  277. }
  278.  
  279. long
  280. _ldiv(a, b) long a, b; {
  281.     int neg = 0;
  282.     long _ldivmod();
  283.     if (a < 0L) { neg++; a = -a; }
  284.     if (b < 0L) { neg++; b = -b; }
  285.     a = _ldivmod(a, b, 0);
  286.     return neg & 1 ? -a : a;
  287. }
  288.  
  289. long
  290. _lmod(a, b) long a, b; {
  291.     long _ldivmod();
  292.     if (a < 0L) a = -a;
  293.     if (b < 0L) b = -b;
  294.     return _ldivmod(a, b, 1);
  295. }
  296.  
  297. long
  298. _ulmul(a, b) unsigned long a, b; {
  299.     long result;
  300.     result = 0L;
  301.     while (a) {
  302.         if (a & 1L)
  303.             result += b;
  304.         a >>= 1L;
  305.         b <<= 1L;
  306.     }
  307.     return result;
  308. }
  309.  
  310. long
  311. _uldiv(a, b) long a, b; {
  312.     long _ldivmod();
  313.     return _ldivmod(a, b, 0);
  314. }
  315.  
  316. long 
  317. _ulmod(a, b) long a, b; {
  318.     long _ldivmod();
  319.     return _ldivmod(a, b, 1);
  320. }
  321.  
  322. long
  323. _ldivmod(a, b, rem) unsigned long a, b; int rem; {
  324.     int i = 0;
  325.     unsigned long result = 0L;
  326.     if (b) {
  327.         while ((a > b) && !(b & 0x80000000L)) {
  328.             i++;
  329.             b <<= 1L;
  330.         }
  331.         while (1) {
  332.             if (a >= b) {
  333.                 a -= b;
  334.                 result++;
  335.             }
  336.             if (i == 0) break;
  337.             i--;
  338.             result <<= 1L;
  339.             b >>= 1L;
  340.         }
  341.         return rem ? a : result;
  342.     }
  343.     else    i/0;    /* should case a "divide zero" trap */
  344. }
  345.